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We report results of a Monte Carlo study of doped, diluted magnetic semiconductors in the low car- 
rier density (insulating) regime. We find that the system undergoes a transition from a paramagnet 
at high temperatures to a ferromagnet at low temperatures. However, in strong contrast to uniform 
systems, disorder effects dominate the entire collective behavior, leading to very unconventional 
properties such as (1) magnetization curves that cannot be described by theories based on critical 
fluctuations, or on spin wave analysis over any significant fraction of the phase diagram; (2) a large 
peak in susceptibilty well below the ordering temperature; and (3) specific heat curves that point 
out the inadequacy of a classical Heisenberg model for spin-5/2 Mn ions. A picture of percolating 
magnetic polarons appears to describe the data well, and leads to a prescription for correcting the 
unphysical results of the classical continuous spin model in terms of a discrete vector model. 
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In diluted magnetic semiconductors (DMS), such as 
Zni_a;Mna;Te and Cdi_xMna;Se, carriers localized at im- 
purities can induce sizable magnetic moments on the 
scale of their hydrogen-like orbits. These entities, known 
as bound magnetic polarons (BMPs), are created by the 
sp-d exchange interactions of localized carriers with mag- 
netic ions in their vicinity. The mechanism and ther- 
modynamics of mutually independent BMPs in lightly 
doped DMS are well understood ||l|,P 
iments on doped, p-type ZnMnTe 
doping densities (but clearly below the transition to a 
metallic state), have shown a significant enhancement in 
susceptibility at low temperatures (T), implying a weak 
ferromagnetic polaron-polaron interaction. Based on an 
idealized, but exactly soluable, model of the single po- 
laron, Wolff et al. El concluded that Mn^"*" spins in the 
intervening region between polarons formed around two 
bound carriers can generate an effective ferromagnetic 
interaction between the two polarons which likely over- 
comes the direct antiferromagnetic interactions between 
impurities (due to virtual carrier hopping, as in usual 
nonmagnetic semiconductors) . 

In this paper, we report a Monte Carlo study of a more 
realistic model (with spatially dependent carrier wave- 
functions) and demonstrate the generation of such ef- 
fective ferromagnetic interactions within a model with 
only antiferromagnetic exchange interactions. We find 
that the system undergoes a transition to a ferromagnetic 
state below T = Tc ] however, the state is characterized 
by a very inhomogeneous spatial magnetization pattern, 
with a large (but finite) peak in the magnetic suscep- 
tibility at a temperature Tp well (one or two decades) 
below Tc, and Tp/Tc can be tuned, depending on the 
carrier density. Furthermore, the classical (continuous 
spin) Heisenberg model is found to yield unphysical re- 
sults for the low-T specific heat, which is corrected by 
using a discrete vector spin model. 

The sp-d exchange interaction of a localized carrier 
with a magnetic ion on a distance R can be described 
by the spin Hamiltonian: 



H = Jo|0(R)ps.S, 



(1) 



where Jq is the exchange constant, while s and S are 
spins of the carrier and the magnetic ion, respectively. 
For this study, we consider for concreteness the case of 
electron doping (s = 1/2); our qualitative results are, 
however, valid for both electron and hole doping. The 
orbital wavefunction of the localized carrier 4){r) is taken 



to be of a hydrogenic form (j){r) = e 



-r/of 



for simplicity 



in which as is the Bohr radius of the hydrogen- like orbit; 
extending to other spatial forms is straightforward. 

On a zincblende lattice, Mn ions are introduced ran- 
domly to replace cations (one of the two interpenetrating 
fee lattices) while the dopants replace anions randomly 
on the other sublattice. The ratio of Mn ions to dopants 
is taken to be 32:1, which corresponds roughly to 0.1% 
Mn at a doping concentration of 10^^cm~^. The Bohr 
radius a^ is taken at a typical value of 20A and the lat- 
tice constant ao = hA. Following Wolff et al. Q, the 
direct coupling between carrier spins has been neglected. 
In fact, we have found that introducing such a coupling 
at the densities studied has no qualitative effects on the 
thermodynamic phase diagram, and very minor quanti- 
tataive changes. We also exclude direct Mn-Mn interac- 
tions among 0.1% substitutional manganese ions, since 
only 1% of these Mn ions have nearest neighbors. (The 
main effect of these pairs or larger clusters is to effec- 
tively reduce the percentage of active Mn spins, since 
most cluster interactions reduce the net magnetic mo- 
ment - e.g. pairs form inert singlets). The Hamiltonian 
can be written as 



n 



i,3 



J(rj,Rj)s., -Sj 



(2) 



where r^ and s^ are position and spin of the i-th carrier 
spin, while Rj and Sj those of the j-th Mn spin. All 
the spins are treated as classical Heisenberg spins. The 
exchange between a Mn spin and a carrier spin is given 
by, 



J(r„R,) = Joe-'l"-*-^^!/'^^. 



(3) 



We carried out simulations on lattices of linear size 
L = 40 to 80, which contain A'c = 8 to 64 dopants and 
Nd = 256 to 2048 Mn ions. We averaged up to 3000 
samples per data point, depending on L and T. Jq is 
taken as the unit of T. The equilibration of each system, 
which consists of as many as 2048 Mn spins, is checked by 
the following technique. We simulate two mutually inde- 
pendent samples with identical locations of particles but 
different initial spin configurations. One sample starts 
from totally random spin configuration, while the other 
from an ordered state, where the electron spins and the 
Mn spins arc oriented fcrromagnetically within each kind 
but opposite to the other kind. The evolutions of the two 
samples correspond to the relaxations from high T and 
low T limit respectively. We expect to obtain the same re- 
sults, such as magnetization and energy, in both samples 
for sufficient long time t, when equilibrium is reached. 

Figure n^ shows the evolution of average magnetization 
(|M|) with t for 256 Mn spins and temperature T — 0.01. 
From the convergence of the two curves we can deter- 
mine the equilibration time for that L and T. Because 
the magnetization values in the two replicas approach the 
equilibrium from different directions, we are able to ob- 
tain the equilibrium results to desired precision without 
wasting precious computer time. 



spins. At large T, g{L) decreases with number of Mn 
spins, or system size, while at low T, it increases with 
system size. The curves seem to cross around T — 0.014, 
at which g{L) is independent of i. We identify the cross- 
ing point as the transition temperature Tc, as suggested 
by finite size scaling theory. 
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FIG. 2. Plot of g{L) for system of 256, 864, and 2048 Mn 
spins at different temperatures around Tc = 0.014. 
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FIG. 1. Plot of magnetization per spin for a system of 256 
Mn spins at T = 0.01 as a function of time. The upper curve 
starts from ferromagnetic spin configuration, while the lower 
one from random configuration. 

To determine the critical temperature in the model, we 
look, in particular, at the following dimensionless quan- 
tity widely used in spin systems |q|: 



9iL) 



. <\M\^> 
'< |AfP>2 



(4) 



for system of various sizes (numbers of Mn and electron 
spins) at different temperatures. The coefficients in the 
above equation have been chosen so that g{L) = 1 at 
T = 0, and as T — > cxd. In Fig. || we plot g{L) as a 
function of temperature T, for different number of Mn 
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FIG. 3. Magnetization per spin versus temperature for 
system of 864 and 2048 Mn spins. 

In Fig. we plot the magnetization per spin at differ- 
ent temperatures for systems of 864 and 2048 Mn spins. 
Plotted on a logarithmic temperature scale, the mag- 
netization curves look similar to what we observe in a 
typical ferromagnet on a linear plot. On a linear plot, 
however, it is seen that (|Af|) does not saturate even 
down to Tc/5. This can be understood as follows. The 
driving force behind the ferromagnetic transition is the 
ordering of Mn spins around each localized carrier into 
magnetic polarons. The radius of the polaron surround- 
ing a dopant ion can therefore be estimated by equaliz- 
ing the exchange energy at such distance to the ambient 



temperature. For localized wavefunctions (in particular 
of the hydrogenic form), therefore, the radius of a po- 
laron increases logarithmically with inverse temperature, 
tt ^ (as/2) ln(JosS'/r). At a certain T, the radius of 
polarons grows large enough so that polarons percolate 
into a macroscopic ferromagnetic cluster, whose forma- 
tion signals the transition, Tc- However, since the per- 
colating cluster does not contain a majority of the spins 
until temperatures which are exponentially lower than 
Tc, most Mn spins remain free even in the ferromagnetic 
phase. It is these uncoupled spins which give rise to the 
slow saturation of magnetization at low T, and to the 
giant susceptibility peak below. The magnetization ap- 
pears to approach the saturation value linearly at low T, 
within the precision of our simulation. 

Figure H shows the plot of magnetic susceptibility as 
a function of T for different system sizes. Susceptibil- 
ity increases dramatically when approaching the critical 
temperature Tc = 0.014, at which it diverges in the ther- 
modynamic limit. Remarkably, these exists a huge peak 
at T = Tp two orders of magnitude below T^,, which ap- 
pears to be L-independent for large enough lattices. Fur- 
ther analysis shows the peak at Tp is associated with 
single spin excitations of those Mn spins far away from 
any of the bound electrons. 



each of these Mn spins contributes S'^ /3T to susceptibil- 
ity in low temperature limit, the total susceptibility can 
be written as: 



Xtotal i^rji ^ 1 



(6) 



which fits the susceptibility peak in the simulation results 
reasonably well as shown in Fig. 0. Since finite system 
size prohibits the existence of manganese spins very far 
from any electron spins, the susceptibility peak at lower 
temperature is suppressed in small systems. Such an ef- 
fect can be seen in the smallest system (256 Mn spins) 
in our simulation. 

In Fig. @, we show the specific heat per spin as a func- 
tion of T. Notably, the specific heat does not vanish at 
low T, approaching unity instead. This is not surpris- 
ing for a classical Heisenberg spin model, since the low 
temperature behavior must obey the classical equiparti- 
tion theorem. In order to mimic the quantum mechan- 
ical characteristics of Mn spins without carrying out a 
full quantum Monte Carlo simulation, we have adopted 
a discrete vector model, in which each Mn spin can be 
oriented only along one of the six [100] directions. Spin 
flips in the vector model cost a finite energy therefore are 
frozen at very low temperatures. 



800 

700 
600 
500 
400 
300 
200 
100 



y ¥ 256 + 

^ xx ^ 864 X 

X * 2048 X 

"" .,, f th 



'".'X -^ 



1e-06 



1e-05 



0.0001 



0.001 

T 



0.01 



0.1 



FIG. 4. Plot of susceptibility versus temperature T for 
systems of various number of Mn spins. Solid curve is the 
estimate of the contributions from manganese spins outside 
the percolating polaron cluster. 

Well below Tc, it is convenient to view the system as 
a cluster of percolating polarons with a large number of 
Mn ions outside the cluster, which are interacting weakly 
with distant electron spins. Modeling the polarons as 
rigid shells of radius vt , the number of Mn ions outside 
the polaron cluster can be estimated: 



Nt - NMP{rT) = A^A/e-^-P^'-T/a^ 



(5) 



where Nm is the total number of Mn ions and pc the 
doping density. P{r) is the probability of a magnetic 
ion having no bound electrons within distance r. Since 
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FIG. 5. Plot of Cv versus temperature in the classical 
Heisenberg spin model. 

Results for a 12-state vector model has been investi- 
gated in more detail and published elsewhere 0]. In the 
two disctrete spin models, the equilibration of a system, 
especially starting with random initial configurations, be- 
comes extremely slow at low temperatures with single 
spin flips, when energy change associated with such a 
flip becomes exponentially larger than the temperature. 
In fact, most of the samples are frozen in metastable 
states. A cluster algorithm, which is capable of updating 
large number of spins (within a polaron or a percolating 
polaron cluster) in a single step, has been introduced to 
speed up the equilibration R. The algorithm has been 



found to attain convergence at low temperatures, espe- 
cially in small samples. 

The introduction of discrete spin models leads to qual- 
itatively similar magnetization and susceptibility curves. 
The magnetic susceptibility shows a peak two orders of 
magnitude below Tc- Yet, the discreteness of energy lev- 
els simulates the quantization of spin energy levels as de- 
sired. Figure |g shows the plot of specific heat of a system 
of 864 Mn spins at different temperatures in the 6- vector 
model. As a result of discrete symmetry, specific heat 
per spin drops to zero in low temperature limit. Follow- 
ing a similar argument as for the susceptibility, one can 
show pi 



Cv ^ pcaBrrpe 



-4Tvpar^/3 



(7) 



which fits the simulation results very well. 
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FIG. 6. Cv at different temperatures for a system of 864 
Mn spins in the discrete vector model. The solid line is the 
theoretical estimate. 

The simulations reported above were for a low con- 
centration of magnetic ions and a low concentration of 
dopants with a ratio such that the interaction between 
Mn spins and between carrier spins can be neglected, 
and further, the dominant carrier-Mn interaction can be 
well estimated through the isolated carrier wavefunction. 
However, as a consequence, the transition temperature is 
very low which makes relevant experimental observations 
almost impossible. We can raise the temperature range of 
interest by increasing the concentration of dopants (rid)', 
however, working within a pure spin model restricts the 
carrier density to lie below the metal-insulator transition 
(MIT) density so the system stays insulating. Increasing 
both the Mn concentration and dopant concentration by 
a factor of 2 gets us close to the MIT {nj as = 0.25). 
Figure shows the susceptibility for system of 512 and 
1728 Mn spins as a function of T. Both T^ and Tp increase 
significantly, still leaving roughly one order of magnitude 
difference between the two temperatures. There is a de- 
crease in the height of the low temperature peak, because 



increasing the dopant concentration reduces the number 
of weakly coupled Mn spins. Our theoretical estimate 
again fits these results very well, which implies that the 
main contribution of the low temperature peak comes 
from such isolated Mn spins. 

While these simulations were done for the insulating 
regime, we expect that the onset of metallicity, especially 
just above the MIT, will not change qualitatively the un- 
usual aspects of the magnetic behavior of doped DMS. 
In fact, unusual magnetization curves have been reported 
in experiments on metallic Gai^xMuxAs, a III-V DMS, 
where Mn contributes both a localized magnetic moment 
and a carrier (but the system appears to be heavily com- 
pensated, so the carrier density is only about 10% of the 
Mn density) j^]. These effects, arising from disorder in 
the Mn and dopant positions, have been left out in recent 
models [^Q that have appeared in the literature since 
the completion of this work. Preliminary results of a nu- 
merical study of metallic doped III-V DMS ||ll[ taking 
into account the disorder in the Mn positions shows that 
qualitative effects seen in this study remain important at 
the densities well beyond the MIT. 
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FIG. 7. Plot of susceptibility versus temperature 

T for systems of various number of Mn spins near the 
metal-insulator transition concentration. Solid curve is the 
estimate of the contributions from Mn spins outside the per- 
colating polaron cluster. 

In summary, we have demonstrated via a Monte 
Carlo study, a ferromagnetic transition in doped, diluted 
magnetic semiconductors in the low doping (insulating) 
regime. Though the model studied has no frustration, 
it has a very high degree of disorder, characterized by a 
wide distribution of exchange couplings. As a result, the 
system exhibits very unusual properties, such as thermo- 
dynamic properties (e.g., magnetization, magnetic sus- 
ceptibility, specific heat) which are very different from 
homogeneous systems. In particular, standard expan- 
sions used for homogeneous systems - critical fluctuations 
around the transition temperature, or spin wave expan- 
sions around the ordered state (T = 0), do very poorly 



in describing the thermodynamics over most of the phase 
diagram. There is considerable thermodynamic entropy 
down to very low temperatures, leading to an unusual 
peak in the magnetic susceptibility and specific heat at a 
temperature well below the transition temperature. Our 
simulations demonstrate the necessity of incorporating 
the discrete nature of the quantum Heisenberg spins (de- 
spite the ferromagnetic nature of the ordered state) , and 
the possibility of unusual thermodynamics in the T ^ 
limit also. Some of the qualitative features are expected 
to survive in the case of metallic doped DMS as well. 
This research was supported by NSF DMR-9809483. 
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